Knee joint biomechanics and cartilage damage prediction during landing: A hybrid MD-FE-musculoskeletal modeling

Understanding the mechanics behind knee joint injuries and providing appropriate treatment is crucial for improving physical function, quality of life, and employability. In this study, we used a hybrid molecular dynamics-finite element-musculoskeletal model to determine the level of loads the knee can withstand when landing from different heights (20, 40, 60 cm), including the height at which cartilage damage occurs. The model was driven by kinematics–kinetics data of asymptomatic subjects at the peak loading instance of drop landing. Our analysis revealed that as landing height increased, the forces on the knee joint also increased, particularly in the vastus muscles and medial gastrocnemius. The patellar tendon experienced more stress than other ligaments, and the medial plateau supported most of the tibial cartilage contact forces and stresses. The load was mostly transmitted through cartilage-cartilage interaction and increased with landing height. The critical height of 126 cm, at which cartilage damage was initiated, was determined by extrapolating the collected data using an iterative approach. Damage initiation and propagation were mainly located in the superficial layers of the tibiofemoral and patellofemoral cartilage. Finally, this study provides valuable insights into the mechanisms of landing-associated cartilage damage and could help limit joint injuries and improve training programs.


1) Model Structure and Simulation:
In the present work, we use a mesoscopic model of the collagen fibril. This model is built based on the geometry of the atomic scale tropocollagen molecules, which represent the constitutive elements of the fibril. The molecule is formed by a triple helix ( two α-1 chains and one α-2 chain), where each chain is composed of a succession of amino acids (glycine, proline, hydroxyproline, hydroxylysine, Arginine, … ). A detailed characterization is found in Gauteri et al., 1 . Furthermore, the availability of molecular structure of tropocollagen 2 and the elementary interactions between atoms forming the molecule (C, N, O, H) 3-6 provide a feasible process to model collagen molecules [7][8][9] via direct (conventional) molecular dynamics (MD) simulations.
Such simulations are also possible because a single molecule is about 300 nm long, 1.5 nm wide, and contains about 42 thousand atoms, which is a reasonable simulation size for MD.
If it is considered to model the whole fibril, the direct molecular approach is no longer an option: In fact, a typical fibril ranges from 20 nm to a few hundred nanometers in diameter and therefore contains 50 million to a few billion atoms. Unfortunately, modeling and simulating such a fibril is not possible using direct MD because of computational limitations. Therefore, the concept of beads or superatoms in modeling collagen molecules was introduced by Buehler et al. 10 . The coarse-grained approach is based on the concept of abbreviating each set of geometrically adjacent atoms into one superatom (or bead) aggregating together to generate a collagen molecule ( fig. 1). These molecules were connected in each proximity via mature (trivalent) and immature (divalent) crosslinks ( fig. 1) to generate the fibril. Thus, this approach significantly reduces the simulation's complexity and allows access to sizes otherwise inaccessible by conventional MD.
In this work, we follow the same bead definition by Buehler et al., 11 , considering 218 beads per molecule (about ~192 atoms per bead) ( fig. 1). Using the exact bead definition is crucial since changing the number of beads will change the characteristics of the interatomic bonds, hence requiring the definition of a new force field. Since the same bead configuration in Buehler et al., 11 was considered, the interatomic potential developed for this specific configuration can be used. Starting from the defined molecule, the fibril is then built by replicating the tropocollagen molecule orthogonally to its principal axis. A hexagonal configuration was used with a distance of 16.52Å between molecules, forming a fibril of 21.5nm diameter and containing 151 molecules    The formulation of the coarse-grained model of collagen molecules was implemented in several studies [11][12][13] and was proven to mimic the actual behavior of the fibril accurately. Three main energies govern the force field.
The interatomic energy ( ) is a pairwise Lennard Jones (LJ) interaction between beads from different molecules responsible for keeping the fibril together in the radial direction. The LJ potential is given by: where and represent respectively the characteristic distance and the minimum energy of the LJ potential, parameters for pairwise potential are given in Table 3. The bond energy is a hyperelastic interaction between two adjacent beads from the same molecule. It is represented by three regimes of potential energy that is given by: where and 1 are spring constants, 1 is the distance at which the hyper-elastic behavior of the bond is triggered, is the bond-breaking distance, and 1 it is a constant calculated to ensure the continuity of the force field. Parameters for the bond energy are given in Table 2. The angular energy is a harmonic three-body interaction between three adjacent super atoms from the same molecule to control the bending angle between the beads.
where Ѳ represents the bending strength, represents the equilibrium angle, and represents the actual angle between the three consecutive beads. Parameters for angle energy are given in Table 3. Enzymatic crosslinks are protein-protein bonds that make up most of the crosslinks in collagen fibril. Enzymatic crosslinks are initially formed between telopeptide and helical residues producing immature (divalent) crosslinks connecting the end of the tropocollagen molecule to the nearest neighbor from an adjacent molecule. Then, this immature crosslink may react with another telopeptide residue producing a mature (trivalent) crosslink joining three collagen molecules by connecting the end of the molecule and the two nearest neighbors from adjacent molecules. In this work, hyperelastic behavior is considered for both divalent and trivalent crosslinks (parameters in Table 2). The ratio of trivalent crosslinks to the total number of enzymatic crosslinks is considered to be 33% 14 . Then, we vary the crosslink content to simulate their effect on the mechanical behavior of the fibril. The coefficient β represents the density of molecule ends connected to beads from other molecules (a coefficient β = 100% corresponds to two connected ends per molecule). The fibril model was created using MATLAB R2021A by averaging the geometric positions of the atoms in the 3HR2 PDB entry and replicating the molecule in the radial directions. All MD simulations were performed using LAMMPS molecular dynamics software 15 .
A 10fs timestep was used. The fibril was relaxed at 300°K for 1 ns using NPT, then NVT to release the residual stress for 1ns. To model the tensile deformation of the fibril, an axial velocity constraint was imposed on both ends of the fibril while keeping the periodic boundary condition in the longitudinal direction to mitigate any surface energy effects. We then use the virial stresses and the fibril volume to compute the stress-strain curves. Finally, the visualization of the results was performed using the OVITO package 16 .

2) Model Output
Figure 3(a) shows the stress-strain curve for a single molecule. The curve can be divided into three zones describing three different mechanisms. Zone I where the molecule was stretched along its principal axis and resulting a small stress that was due to the change of the angle energy only. Then zone II with a uniform stretch of the molecule that led to an elastic constant of ~7.9 GPa. Finally, Zone III with a sharp change in strength occurs as the bond distance reaches the hyper-elastic critical distance r1 (elastic constant is ~47.1GPa). The molecule continues to stretch uniformly until reaching the breaking point. The molecule breaks when interatomic distances reach the breaking distance . Figure 3(b) shows the stress-strain curve for the collagen fibril for different crosslink densities. The overall trend of the curves is consistent with previous work.
Increasing crosslink densities increases both the ultimate tensile strength and the ultimate tensile strain of the fibril. When crosslinks are present in the fibril, they provide additional resistance to the shearing between molecules and therefore retard the sliding threshold, increasing the ultimate tensile strain and stress.

II) Patellar tendon and ligaments model
The hierarchical composite structure of the ligament is illustrated in figure 4. The tissue is considered as a homogenized continuum in the computations presented herein. Each continuum material point reflects the tissue's statistically homogenous representative volume element (RVE) response. As a result, the macroscopic deformation gradient may be defined as the volume average of the deformation gradient over the RVE under homogeneous or periodic boundary conditions.

1) Fibril Model
The first step is to define the multiplicative decomposition of the deformation gradient to illustrate the relationship between the uniaxial and shear deformation 17,18 . Thus, the total deformation gradient tensor = , where s and f stand for shear and uniaxial deformation. This decomposition approach yields a general expression of the strain energy function at the fibril level in the following form: where the shear moduli µ fl is considered as a function of the elastic fibril deformation with The hyperbolic structure of the strain energy function is beneficial in fitting the stiffness evolution of the fibril predicted by MD simulation during molecular analyses [19][20][21][22] . The plastic flow in the fibril was controlled by the effective stress as follow where M is the Mandel stress where and ne stand for the macroscopic Kirchhoff stress and the current fibril direction, respectively. By combing equations (7) and (8), it can be easily shown that effective stress takes the form as follows, The following equation gives the flow resistance 23,24 .
where h, Φ, and Φ are the softening or hardening rate, yield strength, and saturated flow strength of the fibril, respectively. Owing to fibril softening related to the breakage of crosslinks, the saturated flow strength, Φ , is generally chosen to be smaller than the initial yield strength of the fibril (Φ 0 ) 24 . The latter is a function of the crosslink's density ( ) and is defined here by the density function derived from our mesoscopic model ( fig. 6). Solving the equality between the yield strength of the fibril (Φ) and the effective stress, lead to a critical value of 4 (when 4 = cr ) that informing the yielding initiation. After that, the single crystal plasticity model was employed to drive the evolution of the plastic strain in the fibril [23][24][25] : where ̇ is the plastic strain rate, and Φ is the yield strength of the fibril. Here the Karush-Kuhn-Tucker loading/unloading conditions of the system was used as follow After the calculation of the plastic strain rate, the plastic deformation tensor was integrated using the plastic velocity gradient formula 23 Then, using the deformation gradient decomposition equation, the elastic gradient tensor of deformation was updated as follows: Fibrils' plastic stretch was calculated from the basic definition of stretch, = √ 4p (15) where is assumed to be unity throughout the elastic deformation ( = = ). Then, the yield distribution and its beginning were identified when the parameter ( ) become greater than one.

2) Fiber model
The collagen fiber is modeled as fibril reinforced composite material ( fig. 7) with incompressible neo-Hookean matrix characterized by the strain energy function where is the shear modulus of the fiber matrix material. The elastic strain energy of the fiber under extension is given by where 1 = 4 + 2 4 −1/2 , is the fibril volume fraction of the fiber, and is the volume fraction of the matrix material of the fiber. Then the strain energy function characterizing the fiber under shear is given by where 1 = 1 ( ) = ( ) and represent the shear modulus characterizing the effects of shear interactions at the junction of the fibril and matrix [26][27][28] . The total strain energy density of the fiber is therefore written as

3) Tissue model
The former process (mixed formulation, fig. 8) is used to treat the tissue as fiber reinforced composite material ( fig. 9), where the elastic strain energy of the tissue matrix was considered as where is the shear modulus of the tissue matrix material. The elastic strain energy of the tissue under extension is given by where are the matrix and fiber volume fraction, respectively, and = 1 − . Then the strain energy function characterizing the tissue under shear is given by where The total strain energy of tissue is defined by We can further write the strain-energy function of the tissue as    Finally, all materials parameters, multiplicative decomposition, and invariants are listed in the table below.

Materials parameters
Shear Modulus of the Solid matrix ). However, these four parameters were fixed based on our previous calibration process 30 , where a statistical calibration is followed to find the unknown probability distribution function (PDF) ( fig. 11). The fibril parameters were calibrated to strictly map the MD simulation results at different cross links levels using a nonlinear optimization function from Matlab (lsqnonlin) which minimized an objective function f(x) (Eq. 29) in a least-square sense (∑ ( ) 2 )). The optimal objective function compromises the stress predicted via molecular dynamics and the finite element model.
The initial inputs of the eight fibril parameters were selected based on a plausible range of values using the Latin hypercube sampling MATLAB function. Then these initial inputs were iterated via the optimization procedure to ensure an unbiased estimate of the fibril parameters.
Furthermore, as an additional step of verification, the minimum set of outcomes of the optimization process were then used as an input to globalserch solver (gs) to ensure the best fit used to drive the fibril hyperelastic model employed in the subsequent FE simulations. As a outcome, the coefficient of determination R 2 were found equal to 0.938 ± 0.06. Thus, the fitted curves exhibit acceptable fits to the MDS data for the different level of cross links. The different cross link density of the patellar tendon and ligaments were selected based on our previous aggregate probabilistic calibration (Adouni et al.,30 ). In this selection 20% density were assigned to the patellofemoral ligaments, 60% to ACL, PCL, and LCL, and 80% to the MCL and finally 100% to the patellar tendon.

1) Fibril model:
To simulate fibril behavior, a nonlinear constitutive modeling approach developed by Sajjadinia et al., 31 in which the stress of the fibril can be defined as: Where n and are the current direction and logarithmic strain of the fibril, respectively. and are the collagen stiffening coefficients (initial and strain-dependent) and is a depthdepend elastic materials parameter. The collagen networks were defined as primary and secondary bundles (i) of fibrils based on their orientation relative to the articular cartilage depth.
The fibrils were oriented perpendicular to the subchondral junction and turned gradually in the middle zone to become parallel to the articular surface ( fig. 12). The integration of the fibril stress equation with respect to strain in its axial form led to the strain-energy function ( ) 32 2) Cartilage model: The articular cartilage was then modeled using incompressible hyperelastic behavior reinforced by the developed continuum-damage model of the fibril. The Cauchy stress ( ) in the used model was decomposed into a non-fibrillar ( ) and fibrillar ( ) parts as follow: Where F and J are the deformation gradient tensor and the volumetric deformation, respectively. Gm is the shear modulus and is the relative collagen fibril volume fraction. For more details on the formulation of the material, please see prior works 31,43 .  Table 5. The knee was placed in full extension and the scanning process employed a 3D spoiled  The mesh of the model was obtained through a sensitivity analysis, where a maximum of 6% difference in the von-Mises stress was considered (Table. 6).

2) Fibril orientation
For each C3D8R element, a local coordinate system was used to define fibril orientation and implemented by a python script. This script read the connectivity of each element and defined the local cross-sectional plane (x', y') and its normal vector along the local z'-direction ( fig. 14). Figure. 14: Element local coordinate system (x', y', z') employed to define the orientation of the fibril.

3) Interaction and loading analyses
The explicit algorithm was used during all the simulations with short step time (0.01s) to mimic quasi-static analysis. The boundary conditions were applied gradually using an exponential function for the normalized amplitude ( fig. 15) to attain smooth results.  Step time (s)

V) Knee Model Validation and Verification:
A knee axial compressive loading scenario was simulated to test the capability of the  . 19).
The superficial (7.32 MPa) and vertical (6.01 MPa) fibrils were mostly stressed in the structure of the articular cartilage, specifically in the area located under and away from the central loading area, respectively. Low stress has been computed at the random fibril zone. All these computed parameters were in good agreement with the former experimental and modeling investigations [52][53][54][55] ( fig. 16-19).